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We study numerically the Kuramoto-Sivashinsky (KS) equation forced by external white noise in 
two space dimensions, that is a generic model for e.g. surface kinetic roughening in the presence of 
morphological instabilities. Large scale simulations using a pseudospectral numerical scheme allow 
us to retrieve Kardar-Parisi-Zhang (KPZ) scaling as the asymptotic state of the system, as in the 
ID case. However, this is only the case for sufficiently large values of the coupling and/or system 
size, so that previous conclusions on non-KPZ asymptotics are demonstrated as finite size effects. 
Crossover effects are comparatively stronger for the 2D case than for the ID system. 
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The Kuramoto-Sivashinsky (KS) equation is a 
paradigmatic model for chaotic spatially extended 
systems, arising in a variety of physical contexts, 
like thin solid films, interfaces between viscous fluids, 
waves in plasmas and chemical reactions, reaction- 
diffusion systems, or combustion fronts [l|. Actually, 
in its stabilized form, it has been shown to provide 
a generic model for parity-symmetric systems fea- 
turing a bifurcation with a vanishing wave number 
[3]. In the presence of external fluctuations, a natu- 
ral generalization is provided by the noisy KS (nKS) 
equation, that reads 



dh 
'dt 



-iyW^h-ICW^h+-{Wh)^+r]{r,t), (1) 



where r G M.^, we will take t^, /C, A to be positive pa- 
rameters, and ri{r,t) is a Gaussian white noise with 
zero mean and correlations 

{v{r,tUr\t'))=2DS{r-r')S{t-t'). (2) 

Thus, the nKS equation reduces to the determin- 
istic KS (dKS) equation in the D = case. In- 
deed, Eqs. ©-(HI) appear in a wide variety of phys- 
ical contexts, from e.g. step dynamics in epitaxy J3| 
to surface erosion by ion-beam sputtering (IBS) [J], 
or diffusion-limited growth [5|. In these, h{v,t) can 
be thought of as the position at time t of a moving 
front above point r on a reference line or plane, that 
will be the physical image to be used in this work. 
One of the intriguing features of the dKS equation 
is the fact that, at least for d = 1 Q, its large scale 
properties display kinetic roughening in the univer- 
sality class of the (stochastic) Kardar-Parisi-Zhang 
(KPZ) equation, as also occurs for the nKS equa- 
tion |7|, |8| . This links the two seemingly opposite 
phenomena of pattern formation and scale invari- 



ance within the evolution of a single system at ap- 
propriate time and length scales. However, the two- 
dimensional case d = 2 remains controversial: on 
the one hand, there are opposing claims [3, [Ifl] on 
the asymptotics of the dKS equation vs that of the 
2D KPZ equation; on the other hand, the asymp- 
totics of the 2D nKS equation is not well under- 
stood. In d = 1 the nKS equation indeed belongs 
to the KPZ universality class, as borne out from 
numerical simulations [7[ and dynamic renormaliza- 
tion group (DRG) analysis Q. However, as for the 
KPZ equation, the DRG approach is inconclusive 
in d = 2. Numerical results [il| suggest non-KPZ 
asymptotics, contradicting naive expectations based 
on the structure of the RG flow, in which v seems to 
change sign under renormalization as in d = 1 [8| . 

In this work, we revisit the numerical study of 
the 2D nKS equation. Using an improved numerical 
scheme, we perform large scale simulations that al- 
low us to identify KPZ scaling as the asymptotic 
state. However, this only occurs for sufficiently 
large values of the effective coupling in the system 
and/or system size, previous conclusions on non- 
KPZ asymptotics being due to finite size effects. 
Nevertheless, crossovers are comparatively stronger 
for the 2D case than for the ID system, which has 
possibly prevented earlier works from assessing the 
actual hydrodynamic behavior. Our results may 
guide in the assessment of the large scale behavior of 
physical systems described by the 2D nKS equation. 

Note that the nKS system ([I])-© depends on 
a single free parameter; for instance, by rescal- 
ing r -^ (/C/i/)i/2r, t -^ {lC/v'^)t, and h{v,t) ^ 
{D/vy/'^ h{r,t), it can be written as 



dh 
'dt 



= -\7^h-\7'^h + 



V9 



{Vhf+ar,t), (3) 



where g = }?D jv^ and the rescaled noise ^(r, i) has 
zero mean and variance 



(e(r,t)e(r',i')>=2<5(r-r')<5(t-t')- 



(4) 



Thus, we can study the full phase space of the nKS 
equation as a function of only the coupling constant 
g and the lateral system size L. 

Initially motivated by results for IBS, Drotar et al. 
|ll| solved numerically ([I])-© for several values of 
the parameters i/, /C and A. They pointed out the im- 
portance of the coupling g, but were unable to reach 
a well defined asymptotic regime. In fact, they found 
two different scaling regimes, in terms of the expo- 
nents values determined from the behavior of the 
surface roughness W^it) = {iljl?) Y.AK^^ *) - ^?) 
and height-difference correlation function G'(r, t) = 
((/i(r'-|-r,t) — /i(r',i))^); here, bar denotes space av- 
erage and brackets denote ensemble average. Thus, 
in the presence of kinetic roughening |12| , the rough- 
ness scales &sW ^ t^ before reaching the stationary 
state, after which W ^ L"', where /3 and a are called 
growth and (global) roughness exponents, respec- 
tively. Moreover, G{r,t) ~ t'^'^ [G{r,t) - r^"'"-] for 
r ^ i^/^ [r ^ t^^^], where z = a/ 13 is the dynamic 
exponent, and aioc — « for the standard Family- 
Vicsek behavior, while aioc ^ Q^ in the presence of 
so-called anomalous scaling [y|. 

Taking into account the uncertainties of the es- 
timates in [11[ and spurious effects due to the 
crossover between the two scaling regimes found, 
the results obtained in this reference for the early 
stage of the growth process are compatible with the 
Mullins- Herring fixed point [12|. This is consistent 
with the anomalous scaling found for the height- 
height correlation function (see Fig. 11 in [Il| for 
t < 50). For late times, Drotar et al. showed that 
surfaces produced by the nKS equation display scale 
invariance with exponent values for a and /3 in the 
ranges 0.25—0.28 and 0.16—0.21, respectively. These 
values differ substantially from those associated with 
the KPZ universality class for d = 2, even allowing 
for the spread that the latter have [ij]. For the 
sake of homogeneity, we will take as reference values 
for the exponents those we obtain for the 2D KPZ 
equation with the same numerical procedure that 
will be subsequently employed for the nKS equation 
(see below). These are [ll orpz = 0.39 ± 0.01 and 
/3kpz = 0.24 ±0.01. 

An important remark concerns the numerical 
scheme used for the integration of (HI)-©. Drotar et 
al. chose a standard finite-difference discretization 
for space derivatives. Currently it is accepted that 
such a scheme underestimates both the KPZ nonlin- 
earity and the effective coupling g [1^, [l^ . Here we 
opt for a pseudo-spectral scheme that has been suc- 
cessfully used for the numerical integration of local 
[ISl [171 and non-local [Ig stochastic equations fea- 
turing nonlinearities of the KPZ type. Details of this 
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FIG. 1. (Color online) Exponent values of the 2D nKS 
equation as functions of g, in the weak couphng regime. 
Solid black bullets (ordinates on left vertical axis) pro- 
vide values of the roughness exponent a; blue squares 
(ordinates on right vertical axis) provide values of the 
growth exponent /3. The solid black and dashed blue 
lines indicate our reference values for the exponents of 
the 2D KPZ equation, qkpz = 0.39 and /3kpz = 0.24, 
respectively. 



numerical method can be found e.g. in [Tsl. Il9l. |20|| . 

We start by considering parameter values of the 
nKS equation (ll])-© that correspond to relatively 
small coupling values from g = 2-10~^ up to 
g = 2 ■ 10^. We achieve this by tuning A while keep- 
ing other parameters fixed at L — 512, ly = 0.2, 
/C = 2, D = 1, Ax = 2, and Ai = 5 • lO'^, the lat- 
ter being the lattice spacing and the time step, re- 
spectively. For each parameter set, we measure the 
global surface roughness W{t) and the (circular av- 
erage of) the power spectral density (PSD) or height 
structure factor [12[ 5(k) = (ft-k^i-k) as functions of 
time. Here, huit) is the 2D Fourier transform of 
h{x,t) — h{t). In these simulations, observables are 
averaged over 50 noise realizations. The structure 
factor has been shown to feature more clear scaling 
behavior than real-space correlations in the presence 
of crossover effects [2l|, that are expected here. 

By fitting the long time behavior of W{t) prior to 
saturation and the small k = |k| behavior S(k) ^ 



1/A:^"+^ at the stationary state 12|, we find /3 close 
to (log) for g < 10, increasing up to /3 = 0.20±0.01 
for 5 = 2 • 10'^, see Fig. [TJ The roughness exponent 
a is also consistent with (log) for g < 20, after 
which it increases, reaching up to a ~ 0.39 ± 0.01 
for g — 2000, see Fig.[TJ Thus we conclude that, even 
at small couplings, the KPZ nonlinearity is able to 
tame the linear instability in the nKS equation and 
induce kinetic roughening properties. For g ^ 10 
these are in the (2D) Edwards- Wilkinson (EW) uni- 
versality class [l3l, as for the ID nKS case [71, IM]- 
However, for larger coupling values the scaling be- 
havior is neither EW nor KPZ, although the value of 
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FIG. 2. (Color online) Time evolution of the surface 
roughness for different values of the system size L. For 
all simulations we use u = 0.1, /C = 4, A = 20, D — 50, 
Ax = 1.25, and At = 5 ■ 10"^, leading to 5 = 2 ■ lO''. 
Results have been averaged over a different number of 
noise realizations: 15 for L — 2560, 30 for L = 1280, 
and 100 for the remaining values of L in the legend. The 
purple solid line is a guide for the eye, and has slope 
^KPZ = 0.24. 



aioT g — 2000 seems already reminiscent of KPZ be- 
havior. Since EW and KPZ scaling are precisely the 
two meaningful fixed points that are found through 
DRG analysis [g, we believe that the intermediate 
exponent values found in Fig. [T] are to be thought of 
as non- asymptotic behavior due to the finite system 
size of our simulations. This behavior is analogous 
to that of the ID nKS equation 0,lil. 

In order to confirm this interpretation, we need to 
explore larger coupling and/or system size values. 
By increasing g, indeed we have been able to reach 
an asymptotic state in which the critical exponents 
are compatible with those of the 2D KPZ universal- 
ity class. As seen above, the roughness exponent a 
already reaches a KPZ-compatible value already for 
moderate values of g and L. However, /3 approaches 
its asymptotic KPZ value only very slowly. The fact 
that a reaches its asymptotic value earlier (i.e., for 
smaller g values) than /3 has been also reported in 
other studies of crossover phenomena within kinetic 
roughening |2l| , and may be partially accounted for 
by the exact link between the roughness and the 
PSD, W^{L,t) = J2kS{k,t). Unambiguous assess- 
ment of the asymptotic /3 value is only possible for 
large g and L, as seen in Fig. [2l in which W{t) is 
plotted for several system sizes at a fixed large cou- 
pling value g = 2 ■ 10^. This very long crossover 
between pre-asymptotic and asymptotic states hin- 
ders the possibility to reach the strong coupling KPZ 
fixed point for small system sizes and small values 
of the coupling, which applies to the simulations 
previously reported for this system [ll|, see green 
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FIG. 3. (Color online) Qualitative asymptotic scaling 
of the nKS equation (see main text) in [g, L) parameter 
space. In simulations, v = 0.1 and K, — 4, Ax = 1.25, 
and At = 5 • 10""^ are fixed, different values of g be- 
ing achieved by tuning A and D. Red bullets corre- 
spond to pre-asymptotic scaling (type i behavior) while 
blue squares correspond to strong coupling, KPZ scaling 
(type ii and iii behaviors). As examples, points A, B, and 
C are explicitly discussed in [22| . Green diamonds are 
results from [HI ]. The solid line is a reference fit com- 
puted by least squares, separating preasymptotic from 
asymptotic scaling behavior. 



diamonds in Fig. [S] In order to analyze the situ- 
ation in more detail, we have estimated /? with a 
more robust methodology. For kinetic roughening 
systems, we can obtain the roughness exponent and 
the dynamic exponent simultaneously from a stan- 
dard data collapse of the power spectral density [12[ , 
in which we rescalc k — s- kt^^~ and S — s- S'fc^"+^. 
Therefore, an investigation of the strong coupling 
regime has been carried out for representative points 
in the {g,L) plane by collapsing the corresponding 
PSD functions using 2D KPZ exponents, see [22| for 
some specific examples. Note, in our case we are in- 
terested in asymptotic scaling so that collapse with 
such exponents values is only to be expected for the 
smallest fc's in the system. Following this procedure, 
we have been able to identify three different behav- 
iors for a fixed g and increasing substrate size L: i) 
Pre-asymptotic regime in which the collapse of the 
PSD is poor for every value of L in our range (see 
Figs, la-lb in [221, point A in Fig. [3|); ii) regime 
in which only the lowest wave-numbers (very large 
wavelengths) of the PSD collapse with the KPZ ex- 
ponents (Figs. 2a-2b in [22], point B in Fig. [3]); and 
iii) fully developed strong coupling behavior in which 
KPZ asymptotics is reached immediately after the 
exponential growth (due to the linear instability) of 
the surface roughness (Figs. 3a-3b in [22], point C in 
Fig.[3|). Results are qualitatively summarized in Fig. 
131 In this figure we can see that, for small 5 < 10'^ 



and L ~ 1024, we have not yet been able to reach 
the KPZ regime (type i behavior). ActuaUy, in the 
(7 — >■ limit the nKS equation becomes linear (and 
ill-posed) and KPZ behavior does not occur for any 
L value. We then find an intermediate region (for 
10^ < g < 10^ and 512 < L < 1024) in which dy- 
namic crossover behavior occurs between the pre- 
asymptotic regime and KPZ scaling, albeit the lat- 
ter only applies to the largest accesiblc scales (type 
ii behavior). Finally, for g > 10^, 2D KPZ scaling 
is readily observed even for small systems, L ~ 256 
(type iii behavior). 

In summary, our numerical results show that the 
nKS equation is asymptotically in the KPZ univer- 
sality class in two space dimensions, confirming pre- 
vious expectations derived from RG analysis, and 
generalizing known results in d = 1 to one higher 
dimension. This result moreover can guide the in- 
terpretation of large scale experimental and/or nu- 
merical data obtained in the different contexts for 
which this equation appears as a physical model. 
Note that, even in experimental systems, crossover 
effects may hinder observation of actual asymptotic 
behavior at accesible scales, see [l^ and references 
therein. Moreover, and also of practical implica- 
tions, crossover effects are substantially stronger for 
the 2D case than for the ID case, in the sense that, 
fixing all parameter values including the system size 
L, the nKS equation can be already in the KPZ 
asymptotic state for d = 1 while only preasymptotic 
scaling can be measured for d ~ 2. As an example, 
Ueno et al. Q obtain KPZ scaling for the ID nKS 
equation already at g = 20 (and L = 2-10'*), while in 
our 2D case this value of g leads to non-asymptotic 
scaling for any feasible system size, see Fig. [TJ This 
ID vs 2D difference might be due to the particularly 
strong effect that fluctuations have in one dimen- 
sion, which may aid the approach to the stationary 
state for a given parameter set. 

In the context of the controversy on the univer- 



sality class of the deterministic 2D KS (dKS) equa- 
tion, if an effective description of the dKS equation 
by an "equivalent" nKS equation were achieved as 
in the ID case [23, then our results would imply 
that the asymptotic scaling of the 2D dKS equa- 
tion is in the 2D KPZ class. However, such a link 
is not yet available for the d = 2 case [2J|, and in 
the absence of further progress in that direction the 
controversy remains an important open quesion in 
Nonlinear Science. Reflecting on the complexity of 
this problem, one may draw lessons from the case 
of the related Michelson-Sivashinsky (MS) equation, 
that is a model for e.g. flame front propagation |25| . 
Actually, both the KS and the MS equations take a 
very similar shape in k space, thus 



dthk 

where ^"[-1 



A 



= (i.fc"-0™)/ik + -^[(VM1, (5) 

stands for space Fourier transform, 
{n,m) ~ (2,4) for the KS equation and {n,m) = 
(1,2) for the MS equation [l^. The asymptotic 
states of the deterministic MS equation and of the 
noisy MS equation [obtained by adding a noise term 
to the rhs of Eq. ([5]), much like the nKS generalizes 
the dKS equation] are known to be quite different 
j2a |. Nevertheless, this difference turns out to be 
hard to assess in practice, as unavoidable numerical 
noise (round-off errors) in any numerical integration 
of the deterministic MS equation has been seen to 
transform the problem into that of its stochastic gen- 
eralization [2g|. A similar "practical" difficulty in 
telling properties of the deterministic equation apart 
from those of the stochastic generalization may ap- 
ply in the 2D KS context, although whether that is 
the case remains to be seen in the future. 
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